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PREFACE 

The  present  report  by  Miss  Anneli  Leopold, 
prepared  under  Navy  Contract  N6orl-201,  Task 
Order  No.  1,  supplements  Dr.  Friedrichs'  paper 
Formation  and  Decay  of  Shock  Waves  report 
IMM-NYU  158.   It  discusses  the  range  of  validity 
of  the  numerical  approximation  given  in  the 
latter  report  for  the  treatment  of  decaying  shock 
waves. 

Miss  Leopold  was  greatly  assisted  in  the 
numerical  computations  by  Mr.  John  Butler. 

R.  Courant 


DaCAYINQ  SHOCKS 

A  Comparison  of  an  Approximate  Analytic  Solution 
with  a  Finite  Difference  Method. 


INTRODUCTION. 

In  report  IMM-NYU  158,  Formation  and  Decay  of 
Shock  Waves  by  K.  0.  Friedrichs,  an  approximate  method 
was  developed  to  determine  the  behavior  of  a  decaying 
shock  wave  in  a  gas,  i.e.  the  behavior  of  a  shock  wave 
weakened  by  a  rarefaction  wave  overtaking  it.   This 
method  is  expected  to  yield  rather  accurate  results  In 
the  case  of  weak  shocks;  it  is  of  interest  to  deter- 
mine for  what  range  of  shock  strengths  its  use  may  be 
Justified,  and  beyond  what  shock  strength  a  more  exact 
solution  is  required. 

In  order  to  test  the  validity  of  the  "approximate" 
method  by  comparison  with  an  "accurate"  result,  a  corre- 
sponding problem  for  surface  waves  in  shallow  water  is 
treated  first  by  the  "approximate"  method  of  IMM-NYU  158, 
and  then  by  a  finite  difference  scheme  applicable  to 
shocks  of  any  strength.   The  accuracy  of  the  solution 
obtained  by  the  latter  scheme  is  limited  only  by  the  size 
of  the  mesh  used.   The  discussion  of  a  problem  concerning 
water  wave  phenomena  rather  than  gas  dynamics  proper  Is 
advantageous;  for,  while  the  differential  equations  are 
of  the  same  form  in  both  cases,  there  are  no  entropy 
changes  in  water  so  that  the  finite  difference  scheme 
is  simplified  considerably. 
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For  every  shock  strength  considered,  a  comparison 

of  the  results  obtained  by  both  methods  determines  the 

accuracy  of  the  "approximate"  method.   The  shock  strengths 

treated  here.  I.e.  the  ratios  of  height  h^  of  the 

moving  water  to  the  height  of  h^  of  the  still  water 

(see  fig.  1),  are:  1.2,  1.5,  1.8,  3.0  respectively.   The 

results  show  a  very  good  agreement  of  the  two  methods  of 

^1 
solution  for  all  these  cases,  even  for  ^  =  3.0  which  is 

o 

a  stronger  shock  than  those  usually  classed  among  weak 

shocks. 

The  formulation  of  the  physical  problem  Is  as  follows; 

A  shock  wave--or  bore  as  It  Is  also  called--moves  with 

constant  speed  U  Into  still  water  of  uniform  depth  h^. 

The  height  of  the  water  behind  the  shock  Is  h^^,  and  Its 

velocity  Is  u^   (see  fig.  1). 


1. 


a. 


Figure  1 

At  a  certain  time   t  =  0,  a  vertical  wall  Is  suddenly  In- 
serted into  the  moving  water  at  a  point  x  =  0,  so  that 
a  depression  (or  rarefaction)  wave  is  created  which  will 
catch  up  with  the  shock,  since  such  disturbances  are  propa- 
gated with  a  speed  greater  than  that  of  the  shock  (see 
pages  5  and  6).   Thus  the  shock  will  gradually  be  weakened. 


*  Such  a  bore  may  be  created  by  pushing  a  rigid  vertical 
plate  through  the  water  at  constant  speed. 
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Figures  2  and  3  describe  the  physical  situation  at  two 
Instances  after  the  insertion  of  the  wall. 


I 


Plp^ure  2 


UqU 


Figure  3 

In  section  1  of  this  report,  the  "approximate"  method 
will  be  discussed  and  applied  to  the  problem  as  formulated 
here.   In  section  2,  the  finite  difference  method  is  de- 
scribed and  used.   The  graphs  II  and  III  show  the  results 
of  both  methods. 

A  noteworthy  result  of  these  considerations  is  that 
strong  shocks  weaken  very  rapidly  at  first.   This  suggests 
the  following  treatment  of  strong  shocks:   One  uses  the 
finite  difference  method  initially  but  finds,  after  rela- 
tively few  steps,  that  the  shock  has  weakened  sufficiently 
to  justify  the  use  of  the  "approximate"  method  in  con- 
tinuing the  solution.   Thus  the  method  of  IMM-NYU  158  is 
an  extremely  useful  one  even  in  cases  where  shocks  are 
strong. 
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Let  g  be  the  acceleration  due  to  gravity,  h  the 
depth  of  the  water,  x     the  distance  from  the  wall  and  t 
the  time  after  the  Insertion  of  the  wall.   Let  x^  be 
the  distance  between  the  wall  and  the  discontinuity  at 
the  time  the  wall  is  inserted.  I.e.  at  t  =  0.   We,  use 
the  following  definitions  in  order  to  work  with  dimension- 
less  units: 


u  "Vgh   =  particle  velocity. 


XX^       =    X, 


U  Ygh       =  shock  velocity. 


"1^0 

=  ygh     so 

0^ 

=  ^       and 
^o 

\' 

=  t. 

Wr 


Our  initial  conditions   are: 

h(x)   =  h,     for  0  <  X  <  1.      h(x)  =  h^     for  at  >  1 
1  o 

u(x,0)   =  0     for  X  >  1 

The  problem  is  to  find  the  strength,  velocity  and  position 
of  the  shock,  as  functions  of  time. 
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Section  I«  Approximate  method  (applicable  to  weak  shocks). 

If  P  Is  the  density  and  p  the  pressure  of  a  gas  and 
V  the  velocity  relative  to  the  shock,  then  the  mechanical 
shock  conditions  for  a  gas  are: 

<1)  fo%  =  i'l^ 

(2)  Po  -^J=v2  =  p^  +  ^vl 

A  discontinuity  In  water  satisfies  these  same  equations 
if  p  is  interpreted  as  dgh  (where  d  is  a  constant  de- 
noting the  density  of  water),  and  if  p  is  interpreted 
dh  .   Then  relation 

(»'  ^  =  231^" 

holds.   It  is  the  analogue  of  the  adiabatic  equation  of 
state.   Relation  (3)  does  not  Involve  entropy,  and  the 
adiabatic  exponent  is  2. 

From  these  relations,  the  following  equations  are 
derived: 

(4)  u  -  U  =  -  ^  (1  +  c^) 

(5)  U^  =  -I  c^d  +  c^) 


(6)         u  =  ^^  (c^-  1)  Vl  +  c^  =  4  (c^  -  1) 

cY2"  c 

Equations  (5)  and  (6)  determine  the  shock  velocity  U 

and  the  particle  velocity  u  on  the  left  of  the  shock 

in  terms  of  c  =1/—   on  the  left  side  of  the  shock. 

'^o 
It  was  mentioned  previously  that  a  rarefaction  wave 

was  created  by  the  insertion  of  the  wall.   Since  the  front 

of  a  rarefaction  wave  moves  with  sound  speed  c,   relative 
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to  the  medium  in  front  (whose  velocity  is  u^),  its  equa- 
tion is  -f  =  Ut  +  c-   and  it  catches  up  with  the  shock 
whose  velocity  is  subsonic  relative  to  the  medium  behind 
it.   The  point  of  intersection  x  =  ^^ ,  t  =  k  of  the  shock 


line     U  = 


^-1 


with  the   front  of   the  wave     "^  "  ^1  "*"  *^i 


is 


determined: 
(7) 

(8) 


/    = 


k  = 


^1  -^  ^ 


\ 


> 


u^+c  -U 


/ 


Prom  now  on,  the  point  '^,k  will  be  the  new  origin  of  the 
new  X,  t-coordlnate  system  given  by  x  =  x+^,t  =  t  +  k, 
see  figure  4. 


►X 


Figure  4 


Formulae  (7)  and  (8)  introduce  a  slight  inconsistancy. 
In  the  approximate  method,  the  Riemann  invariant  is 
assumed  to  hold  across  the  shock,  i.e.   u  -  2c  =  -2 
everywhere.   However,   u,   and  c,   were  calculated  di- 
rectly and  u-  -  2e,   is  slightly  different  from  -2. 
The  point  Z   ,k  was  obtained  from  this  u,   and  c^ ,  hence 

differs  slightly  from  the  point  X    ,k   which  would  be 
obtained  by  using  u,  -  Zc-^   =  -2.   However,  this  disagree- 
ment is  small  even  in  case  of  rather  strong  shocks.   The 
two  methods  of  finding  u,   and  c,   enable  us  to  deter- 
mine how  accurate  the  approximate  method  is  Initially. 
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It  Is  one  of  the  essential  features  of  this  approxi- 
mate theory  that  the  reflection  of  the  rarefaction  wave 
from  the  shock  Is  Ignored  (since  Its  effect  enters  only 
In  terms  of  higher  than  second  order).   In  the  x, t -plane 
the  rarefaction  wave  Is  a  centered  simple  wave  (see  Manual 
on  Supersonic  Flow  and  Shock  Waves,  AMP  Report  38. 2R, 
Chapter  III,  Article  21).   As  such.  It  has  the  property 
that  one  family  of  characteristics  consists  of  straight 
lines  through  the  point  x  =  t  =  0  along  which  u  and 
c  are  constant.   Each  of  these  straight  lines  Intersects 
the  negative  x-axls  at  a  point  x  =  '  . §  To  each  value  of 
this  parameter  $  ,  there  belongs  one  and  only  one  straight 
characteristic;  hence,   u  and  c  may  be  considered  func- 
tions of  5  .   Specializing  the  relation  -  =  u  +  c  to 
t  =  k  and  x  =  h  +Z,    one  obtains  2^  =  u(  ^  )  +  c(  5  ). 
On  the  other  hand,  using  the  Riemann  Invariant  across  the 
rarefaction  wave,  one  has  u( §  )  -  2c( 5  )  =  u^  -  2c^. 
Solving  these  two  equations  simultaneously,  one  obtains 
the  values  of  u  and  c  as  functions  of  5  s 

(9)  u(n  =§.[2  i^  +  (u^-2c^)] 

(10)  c(^  )  =  J  [   ^  -  (u^-2c^)] 

If  the   assumption  that   the   Riemann  invariant     u   -  2c 
is   constant   across   the  shock   (see  IMM-NYU   158)   is  used  , 
then  equations    (9)    and   (10)   become: 

(9*)  u(5  )    =§  [2  i^  -  21 

do')  c(!>  )  =i  [    i^  +  2] 


*  This  assumption  is  consistent  with  the  "approximate" 
theory  of  report  IMM-NYU  158,  because  the  change  In  the 
Riemann  invariant  across  the  shock  appears  only  in  third 
and  higher  order  terms. 
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In  order  to  find  the  shock  curve  U  =  f[x(  5  ),'t(l  )], 
a  differential  equation  for  T( 5  )  is  set  up  and  solved. 
We  have  now,  on  the  left  side  of  the  shock, 

i-=-i-  =  u(  S)  +  c(M  or   X  =  [u(  !>)  +  c(S-)]  ^  +  §  ; 


differentiating. 


dx 


dt  .  -rr.  ' 


[u(  5  )  +  c(  S  )]  ^  +  t[u  (  5 )  +  c  (  5  ))  +  1  , 


where  the  prime  denotes  differentiation  with  respect  to  5  • 


But 


dx  _  dx  dt   n  ^^ 
d5  "  dt  dj  -  "  dl 


u  II  =hx{l  )}  ||  +  "t  [u  (5)  +  c  (S)]  +  1 


(11)   or  ^[y.(^(l)^ci$))]=:-tU\l)^c\l))-^l 


n 


Equation  (S()  is  an  ordinary  differential  equation  in  t(  1  ). 


Setting    r   ^4^Lll£-^7Ti  *i3  =  ^  (3>  ),  the  soluti 

/  u-[u(S)+c(5)]  -^ 

I/O 


on 


(12) 


(13) 


t(  S  )  =  e 


-  i(M 


e-I(5) 
U-.(u+c)  ^^5 


\ 


> 


x( 5  )  =  5  +  (u  +  c)t 


/ 


yields  a  pair  of  parametric  equations  for  the  shock  curve. 

In  order  to  obtain  a  simple  explicit  expression  for 

"t(  S),  an  expansion  of  the  quantities  u  and  U  in  terms 

of  the  parameter  cr(5)  =u(!>)  +  c(S)  -  1  (similar  to 


-   9    - 


that  described  In  IMM-NYTJ   158)    Is   made;    terms  up   to   second 
order  are  kept   and   equations    (12)    and    (13)   become: 


(14) 


(15) 


1(5)    =   k 


ff^(36/5   -  ef 
_  o''^(36/5  -  a^)"^ 


-1 


:(5  )    =     5    +    (OT   +   1)    t(  !)) 


where     S     =    5(0)   =  u(0)   +  c(0)    -   1  =  u^    +   c.    -  1   . 

O  XI 

In  terms  of  the  parameter  5  ,  these  equations  can  be  written 
In  the  form        ^ 

-2   1 


r 


(16) 


t(  !>)  =  k 


< 


(41/5 


(41/5  -  8)(^ 

x(5)  =  3  +  -^  t(!,) 


-  •^4^)(3  -  1) 

K^  - 1) 


-  1 


where  s  =  u.  +  c  |  or.  In  terms  of  the  original  coordinates 
X  and  t. 


t(3)  =  k 


(41/5  -  -^-^)(s  -  1) 
(41/5  -  s)  (jL^  -  1) 


(17) 


x(5  )  =  ^^  t(3) 


which  is  convenient  for  numerical  work.   Equations  (17) 
were  used  to  compute  the  shock  curves  represented  on  the 
graph  by  broken  curves. 

It  Is  interesting  that  only  a  part  of  the  depression 
wave  interacts  with  the  shock;  the  remaining  part  leaves 
the  shock  unaffected.   In  the  physical  situation,  the  active 
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portion  of  the  depression  wave  Is  the  forward  moving  part 
whose  height  Is  greater  than  h  (see  figure  5).  In  the 
X, t-plane,  the  active  region  Is  that  In  which  the  straight 
characteristics  satisfy  the  Inequality  -^  =  u  +  c  >  1 
(see  figure  6).   This  behavior  may  be  deduced  by  examining 
equations  (17).   When  ^^'^   =  1,  x  and  t  (the  coordinates 
of  the  point  of  Intersection  of  the  shock  with  the  charac- 
teristic where  5  =  J>  =  k  -  ^)  become  Infinite.   Insert- 


ing this  value  of 


i+J 


In  (9  )  and  (10  ),  one  has: 


u(  3*)  =  0,   c(  5*)  :*  1 


This  illustrates  that  only  the  forward  moving  part  of  the 
water  affects  the  shock  (see  figures  5  and  6),  and  that 
it  alone  decays  the  shock  for  an  infinltaly  long  time. 


Figure  5 


•►X 


Figure  6 
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Suppose  that.  Instead  of  moving  the  wall  with  such  a 
velocity  that  cavitation  occurs  (i.e.  the  case  in  which 
c  varies  from  c.   to  0  across  the  rarefaction  wave), 
the  wall  Is  given  a  constant  velocity  w  greater  than 
u.  -  2c.   hut  less  than  u,.   Then  the  rarefaction  wave 
Is  narrowed  and  a  constant  state  appears  to  the  left  of 
the  variable  zone.   The  particle  velocity  in  this  con- 
stant state  is  just  the  wall  velocity  w.   If  w  <  0, 
nothing  essentially  different  happens,  because  then,  the 
rarefaction  wave  is  still  wide  enough  to  Include  the 
"active"  portion  which  decays  the  shock  for  an  infinitely 
long  time.  If  however,  0  <  w  <  u.,  there  is  an  interaction 
which  ends  after  a  finite  time  during  which  the  shock  de- 
cays from  one  constant  shock  to  a  weaker  constant  shock, 
(see  figures  7  and  8).   The  final  constant  shock  is,  of 
course,  constant  only  within  the  approximation  considered 
here. 


UalL 


Bzfoitdsceiy 


h, 


I 


Uoa 


L 


T 
i_ 


h. 


Figure  7 
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Figure  8 


In  order  to  Investigate  at  which  time  a  shock  has 

decreased  from  h.   to  an  arbitrary  height 

2        ^ 
h  =  c  ,(h  <  h  <  h^),  one  uses  the  equation 


t  [S  (c*)]  =  k 


(17/5  -  e")(s  -  1) 
(41/5  -  s)(c*  -  1) 

.2 


for  t  as  a  function  of  c  =  h  • 

(t  =  k  represents  the  time  at  which  the  decay  begins )> 

Graph  I  represents  the  height  as  a  function  of 
time  in  case  of  the  following  initial  shock  strengths: 

2   ^1 
c!'  =  Tr=  =  1.2  ,  1.5  ,  1.8  ,  3.0  ,  respectively. 

Q 


It  Is  apparent  that  the  height  decreases  rapidly  at  first 
and  then  approaches  Its  asympotlc  value  slowly. 


7o     t 
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Sactlon  II.   Finite  difference  method. 


The  differential  equations  describing  the  flow  In 
this  problem  are: 


N 


(1) 


u.  +  uu  +  2cc  =  0 

t  X  X 


2c.  +  cu  +  2uc  =  0 
t  X  X 


The  corresponding  characteristic  equations  are 


(2)      ^  =  u*o 


u  +  2c  =  const. 


\ 


along  C   characteristics 


(3) 


11  =  '"-  = 


u  -  2c  =  const 


along   C   characteristics 


Replacing  the  derivative  ^  In  (2)  by  the  difference 
quotient,  one  obtains,  for  any  Interior  point  of  the  net 
shown  In  figure  9: 

J^~- )s.s 


Figure  9 
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iA\  ^Ik    "    ^1-1,  k    _    1     r„  ,,  , 


-    X 


(5)  Ik  1,^-1  =  I  [T        +   T  ] 

hk       \,k-l        '^        ^^  ^*^  ^ 

where     S^^^  =  ^x^^  +  c^^^     and     T^^^  =  u^^^  -  c^^^ 
The  slBiultaneoua    solution  of    (4)    and   (5)   yields: 

^,       ^<^i,k-r^i-i,k^^V-i.k^^ik^^-i.k^-^,k..i^^ik-^^i,k-i) 
^^^  ^^k'^i-i,k^-^^ik*^i,k-i' 


(6)< 


^Ik  =  ^i-l,k  *  S    ^*ik   -    *i-l,k^    ^^ik  ■"  ^l-l,k\ 


S     and     T     may  be  expressed  in  terms   of     u,,    ,   u  ,    , 
c^.      and     c..     (by  using   equations    (3))   in  the  follow! 

way: 

f^ik  =  ^k  -^  ^k  =  I    ^^  *  ^^hk  -^  ?   <^   -  2c)ii 


(vK 


^"-  =  ^k-^k=l   (u  +  2c)3_j^^|(u   -20)^^ 


\ 


^ik 


Since  u,.   and  c,.   are  known,  all  interior  points  on  the 
i    line  could  be  found,  provided  u, -   and  c    were  also 
known.   These  together  with  U  .   however,  can  be  obtained 
from  the  relations: 


Solution  of  equations  (8),  (9),  and  (10): 

2  I 1 

A  value   of  U     is   chosen.      Then  c     =     "1/211  +  1/4   -  -x     and     c 

can  be   found.      Next     u     is   determined  from     u  =  K  -  2c, 
where     K  =  u..    +  2c,,    .      Finally     U     is   computed  by 

MO 

U  =  =1 ^  using  the  u  and  c   just  found.   Then  the 

c'^-l 

chosen  U  and  the  computed  U  are  compared  in  order  to 

determine  the  U  for  the  next  step  in  the  iteration. 
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(8)      u.-  +  2Cj^j^  "  ^1  1  ■*"  ^°1  1      ^^°^   ^^^ 


(9) 


u 


il 


°il   *  ^ 


^1 


=  0 


(10) 


u 


11 


^1  ^1 

2    -, 
^1   -  ^ 


(10  )    u 


11 


^1  -^  1 


'11 


(9),  (10),  and  (11)  are  obtained  from  the  shock  conditions, 
(see  page  5). 

Graphs  II  and  III  represent  the  Interaction  of  the 
shock  and  the  rarefaction  wave  In  the  x,t-plane.   The 
solid  curves  and  their  Intersection  points  were  obtained  by 
the  finite  difference  method  just  described.   The  broken 
curve  represents  the  path  of  the  shock  as  obtained  from 
the  approximate  method.   For  x  <  7,  the  shock  lines  ob- 
tained by  the  two  methods  coincide  within  the  accuracy  of 
the  graph.   In  the  first  case,  a  ratio  of  heights  of  1.8 
was  assumed,  while  in  the  second,  it  was  taken  as  3.0  . 
Although  these  ratios  are  much  higher  than  those  usually 
considered  in  weak  shocks  the  agreement  between  the  approxi- 
mate -  and  finite  difference  methods  is  extremely  good. 

For  large  values  of  x  and   t,  the  shock  curves  no 
longer  coincide.   However,  since  the  mesh  of  the  finite 
difference  scheme  is  very  large  in  that  region,  this  de- 
viation may  be  due  to  the  Inaccuracy  in  this  scheme, 
particularly  since  the  other  method  should  gain  in  accuracy 
as  the  shock  becomes  weaker. 
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One  may  conclude  that  even  higher  than  second  order 
terms  In  the  expansion  of  quantities  across  the  shock 
differ  little  from  what  they  would  be  in  the  case  of  a 
continuous  wave.   It  therefore  seems  reasonable  to 
apply  the  method  of  IMM-NYU  158  even  in  cases  where  the 
pressure  ratios  are  rather  high. 
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